####################################################################
## author:    Michael L. Wicki
## contact:   michael.wicki@istp.ethz.ch, ETH Zurich
## file name: mobility_conjoint.R
## Context:   Data Analysis
## started:   2018-07-26
####################################################################

rm(list = ls())
getwd()
setwd("\\\\gess-fs.d.ethz.ch\\home$\\wimi\\Documents\\Sufficiency\\Conjoint")
set.seed(42)
.libPaths( "\\\\gess-fs.d.ethz.ch\\home$\\wimi\\Documents\\R\\win-library\\3.3")

################################ load packages ################################ 

#install.packages("Matrix")
#install.packages("ggthemes")
#install.packages("cjoint")
#install.packages("ggplot2")
#install.packages("effects")
#install.packages("lme4")
#install.packages("data.table")
#install.packages("dplyr")
#install.packages("nFactors")
#install.packages("texreg")
#install.packages("tidyr")
#install.packages("forcats")
#install.packages("Rmisc")
#install.packages("cowplot")
#install.packages("scales")
#install.packages("cregg")



library(cjoint)


################################ load dataset & design ################################ 

load("conjoint_design.RData")

###################### summary statistics sample ##############
sampledescr <- df_conj[c("country", "age", "age_grp", "gender", "educ", "incgrp")]
table(sampledescr$educ,sampledescr$country)
table(sampledescr$incgrp,sampledescr$country)
table(sampledescr$age_grp,sampledescr$country)
table(sampledescr$gender,sampledescr$country)


######################### Conjoint Choice All ######################### 
c_base <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
                 attrib4_lab + attrib6_lab, data = df_conj, 
               cluster=TRUE, respondent.id="id", 
               design=design, baselines = baselines_base)

plot(c_base)

c_attribute1 <- c(0, c_base$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base <- data.frame(c,s)
df_c_base$att <- c(c_base$attributes$attrib1lab[1], c_base$attributes$attrib1lab[2:3], 
                   c_base$attributes$attrib2lab[1], c_base$attributes$attrib2lab[2:5],
                   c_base$attributes$attrib7lab[3], c_base$attributes$attrib7lab[2], c_base$attributes$attrib7lab[1],
                   c_base$attributes$attrib3lab[5], c_base$attributes$attrib3lab[4], c_base$attributes$attrib3lab[3], c_base$attributes$attrib3lab[2], c_base$attributes$attrib3lab[1], 
                   c_base$attributes$attrib5lab[3], c_base$attributes$attrib5lab[2], c_base$attributes$attrib5lab[1],
                   c_base$attributes$attrib4lab[3], c_base$attributes$attrib4lab[2], c_base$attributes$attrib4lab[1], 
                   c_base$attributes$attrib6lab[3], c_base$attributes$attrib6lab[2], c_base$attributes$attrib6lab[1])
df_c_base$att <- factor(df_c_base$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                        "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                        "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                        "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                        "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                        "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                        "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                        labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                         "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                         "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                         "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                         "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                         "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                         "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base$poltyp <- factor(ifelse(df_c_base$att == "Baseline: No new tax" |  df_c_base$att == "Increasing prices by 15%" |  df_c_base$att == "Increasing prices by 30%" | 
                                    df_c_base$att == "Baseline: General government budget" | df_c_base$att == "Public environmental and climate protection programs" | df_c_base$att == "Public programs for low-income households" | df_c_base$att == "Reduce income taxes" | df_c_base$att == "No tax revenues" |
                                    df_c_base$att == "Baseline: Keeping subsidies at current level" | df_c_base$att == "Halving subsidies" | df_c_base$att == "Eliminating subsidies",
                                  "Market-based / push", ifelse(
                                    df_c_base$att == "Baseline: Never restricted" | df_c_base$att == "1 day per week" | df_c_base$att == "3 days per week" | df_c_base$att == "5 days per week" | df_c_base$att == "Always restricted" | 
                                      df_c_base$att == "Baseline: No new requirements" | df_c_base$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                    "Non-market-based / push", ifelse(
                                      df_c_base$att == "Baseline: No campaigns" | df_c_base$att == "Occasional campaigns" | df_c_base$att == "Frequent campaigns", 
                                      "Non-market-based / pull", ifelse(
                                        df_c_base$att == "Baseline: No increase of support" | df_c_base$att == "Reducing prices by 15%" | df_c_base$att == "Reducing prices by 30%",
                                        "Market-based / pull", "Other Policy Features")))))

df_c_base$policy <- ifelse(df_c_base$att == "Baseline: No new tax" |  df_c_base$att == "Increasing prices by 15%" |  df_c_base$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base$att == "Baseline: General government budget" | df_c_base$att == "Public environmental and climate protection programs" | df_c_base$att == "Public programs for low-income households" | df_c_base$att == "Reduce income taxes" | df_c_base$att == "No tax revenues", "2", ifelse(
    df_c_base$att == "Baseline: Keeping subsidies at current level" | df_c_base$att == "Halving subsidies" | df_c_base$att == "Eliminating subsidies", "3", ifelse(
      df_c_base$att == "Baseline: Never restricted" | df_c_base$att == "1 day per week" | df_c_base$att == "3 days per week" | df_c_base$att == "5 days per week" | df_c_base$att == "Always restricted", "4", ifelse(
        df_c_base$att == "Baseline: No new requirements" | df_c_base$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base$att == "Baseline: No campaigns" | df_c_base$att == "Occasional campaigns" | df_c_base$att == "Frequent campaigns", "6", ifelse(
            df_c_base$att == "Baseline: No increase of support" | df_c_base$att == "Reducing prices by 15%" | df_c_base$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base$policy = factor(df_c_base$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))



p_c_base <- ggplot(df_c_base, aes(x=att, y = c, shape = poltyp)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_c_base


#########################  Conjoint Rating All ######################### 
r_base <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
                 attrib4_lab + attrib6_lab, data = df_conj, 
               cluster=TRUE, respondent.id="id", 
               design=design, baselines = baselines_base)

plot(r_base)

c_attribute1 <- c(0, r_base$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base <- data.frame(c,s)
df_r_base$att <- c(r_base$attributes$attrib1lab[1], r_base$attributes$attrib1lab[2:3], 
                   r_base$attributes$attrib2lab[1], r_base$attributes$attrib2lab[2:5],
                   r_base$attributes$attrib7lab[3], r_base$attributes$attrib7lab[2], r_base$attributes$attrib7lab[1],
                   r_base$attributes$attrib3lab[5], r_base$attributes$attrib3lab[4], r_base$attributes$attrib3lab[3], r_base$attributes$attrib3lab[2], r_base$attributes$attrib3lab[1], 
                   r_base$attributes$attrib5lab[3], r_base$attributes$attrib5lab[2], r_base$attributes$attrib5lab[1],
                   r_base$attributes$attrib4lab[3], r_base$attributes$attrib4lab[2], r_base$attributes$attrib4lab[1], 
                   r_base$attributes$attrib6lab[3], r_base$attributes$attrib6lab[2], r_base$attributes$attrib6lab[1])
df_r_base$att <- factor(df_r_base$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                        "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                        "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                        "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                        "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                        "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                        "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                        labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                         "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                         "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                         "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                         "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                         "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                         "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base$poltyp <- factor(ifelse(df_r_base$att == "Baseline: No new tax" |  df_r_base$att == "Increasing prices by 15%" |  df_r_base$att == "Increasing prices by 30%" | 
                                    df_r_base$att == "Baseline: General government budget" | df_r_base$att == "Public environmental and climate protection programs" | df_r_base$att == "Public programs for low-income households" | df_r_base$att == "Reduce income taxes" | df_r_base$att == "No tax revenues" |
                                    df_r_base$att == "Baseline: Keeping subsidies at current level" | df_r_base$att == "Halving subsidies" | df_r_base$att == "Eliminating subsidies",
                                  "Market-based / push", ifelse(
                                    df_r_base$att == "Baseline: Never restricted" | df_r_base$att == "1 day per week" | df_r_base$att == "3 days per week" | df_r_base$att == "5 days per week" | df_r_base$att == "Always restricted" | 
                                      df_r_base$att == "Baseline: No new requirements" | df_r_base$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                    "Non-market-based / push", ifelse(
                                      df_r_base$att == "Baseline: No campaigns" | df_r_base$att == "Occasional campaigns" | df_r_base$att == "Frequent campaigns", 
                                      "Non-market-based / pull", ifelse(
                                        df_r_base$att == "Baseline: No increase of support" | df_r_base$att == "Reducing prices by 15%" | df_r_base$att == "Reducing prices by 30%",
                                        "Market-based / pull", "Other Policy Features")))))

df_r_base$policy <- ifelse(df_r_base$att == "Baseline: No new tax" |  df_r_base$att == "Increasing prices by 15%" |  df_r_base$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base$att == "Baseline: General government budget" | df_r_base$att == "Public environmental and climate protection programs" | df_r_base$att == "Public programs for low-income households" | df_r_base$att == "Reduce income taxes" | df_r_base$att == "No tax revenues", "2", ifelse(
    df_r_base$att == "Baseline: Keeping subsidies at current level" | df_r_base$att == "Halving subsidies" | df_r_base$att == "Eliminating subsidies", "3", ifelse(
      df_r_base$att == "Baseline: Never restricted" | df_r_base$att == "1 day per week" | df_r_base$att == "3 days per week" | df_r_base$att == "5 days per week" | df_r_base$att == "Always restricted", "4", ifelse(
        df_r_base$att == "Baseline: No new requirements" | df_r_base$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base$att == "Baseline: No campaigns" | df_r_base$att == "Occasional campaigns" | df_r_base$att == "Frequent campaigns", "6", ifelse(
            df_r_base$att == "Baseline: No increase of support" | df_r_base$att == "Reducing prices by 15%" | df_r_base$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base$policy = factor(df_r_base$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_r_base <- ggplot(df_r_base, aes(x=att, y = c, shape = poltyp)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_r_base


#########################  Conjoint Rating Controls #########################
r_base_ctrl <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                      attrib5_lab +  attrib4_lab +  attrib6_lab + age + country + gender + educ + stopdriving, data = df_conj, 
                    cluster=TRUE, respondent.id="id", 
                    design=design, baselines = baselines_base)

plot(r_base_ctrl)

c_attribute1 <- c(0, r_base_ctrl$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_ctrl$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_ctrl$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_ctrl$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_ctrl$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_ctrl$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_ctrl$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_ctrl$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_ctrl$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_ctrl$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_ctrl$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_ctrl$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_ctrl$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_ctrl$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_ctrl$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_ctrl$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_ctrl$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_ctrl$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_ctrl$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_ctrl$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_ctrl$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_ctrl$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_ctrl$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_ctrl$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_ctrl$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_ctrl$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_ctrl$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_ctrl$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_ctrl$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_ctrl$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_ctrl$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_ctrl$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_ctrl$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_ctrl$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_ctrl$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_ctrl$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_ctrl <- data.frame(c,s)
df_r_base_ctrl$att <- c(r_base_ctrl$attributes$attrib1lab[1], r_base_ctrl$attributes$attrib1lab[2:3], 
                        r_base_ctrl$attributes$attrib2lab[1], r_base_ctrl$attributes$attrib2lab[2:5],
                        r_base_ctrl$attributes$attrib7lab[3], r_base_ctrl$attributes$attrib7lab[2], r_base_ctrl$attributes$attrib7lab[1],
                        r_base_ctrl$attributes$attrib3lab[5], r_base_ctrl$attributes$attrib3lab[4], r_base_ctrl$attributes$attrib3lab[3], r_base_ctrl$attributes$attrib3lab[2], r_base_ctrl$attributes$attrib3lab[1], 
                        r_base_ctrl$attributes$attrib5lab[3], r_base_ctrl$attributes$attrib5lab[2], r_base_ctrl$attributes$attrib5lab[1],
                        r_base_ctrl$attributes$attrib4lab[3], r_base_ctrl$attributes$attrib4lab[2], r_base_ctrl$attributes$attrib4lab[1], 
                        r_base_ctrl$attributes$attrib6lab[3], r_base_ctrl$attributes$attrib6lab[2], r_base_ctrl$attributes$attrib6lab[1])
df_r_base_ctrl$att <- factor(df_r_base_ctrl$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                                  "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                                  "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                                  "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                                  "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                                  "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                                  "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                             labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                              "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                              "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                              "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                              "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                              "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                              "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_ctrl$poltyp <- factor(ifelse(df_r_base_ctrl$att == "Baseline: No new tax" |  df_r_base_ctrl$att == "Increasing prices by 15%" |  df_r_base_ctrl$att == "Increasing prices by 30%" | 
                                         df_r_base_ctrl$att == "Baseline: General government budget" | df_r_base_ctrl$att == "Public environmental and climate protection programs" | df_r_base_ctrl$att == "Public programs for low-income households" | df_r_base_ctrl$att == "Reduce income taxes" | df_r_base_ctrl$att == "No tax revenues" |
                                         df_r_base_ctrl$att == "Baseline: Keeping subsidies at current level" | df_r_base_ctrl$att == "Halving subsidies" | df_r_base_ctrl$att == "Eliminating subsidies",
                                       "Market-based / push", ifelse(
                                         df_r_base_ctrl$att == "Baseline: Never restricted" | df_r_base_ctrl$att == "1 day per week" | df_r_base_ctrl$att == "3 days per week" | df_r_base_ctrl$att == "5 days per week" | df_r_base_ctrl$att == "Always restricted" | 
                                           df_r_base_ctrl$att == "Baseline: No new requirements" | df_r_base_ctrl$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_ctrl$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                         "Non-market-based / push", ifelse(
                                           df_r_base_ctrl$att == "Baseline: No campaigns" | df_r_base_ctrl$att == "Occasional campaigns" | df_r_base_ctrl$att == "Frequent campaigns", 
                                           "Non-market-based / pull", ifelse(
                                             df_r_base_ctrl$att == "Baseline: No increase of support" | df_r_base_ctrl$att == "Reducing prices by 15%" | df_r_base_ctrl$att == "Reducing prices by 30%",
                                             "Market-based / pull", "Other Policy Features")))))

df_r_base_ctrl$policy <- ifelse(df_r_base_ctrl$att == "Baseline: No new tax" |  df_r_base_ctrl$att == "Increasing prices by 15%" |  df_r_base_ctrl$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_ctrl$att == "Baseline: General government budget" | df_r_base_ctrl$att == "Public environmental and climate protection programs" | df_r_base_ctrl$att == "Public programs for low-income households" | df_r_base_ctrl$att == "Reduce income taxes" | df_r_base_ctrl$att == "No tax revenues", "2", ifelse(
    df_r_base_ctrl$att == "Baseline: Keeping subsidies at current level" | df_r_base_ctrl$att == "Halving subsidies" | df_r_base_ctrl$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_ctrl$att == "Baseline: Never restricted" | df_r_base_ctrl$att == "1 day per week" | df_r_base_ctrl$att == "3 days per week" | df_r_base_ctrl$att == "5 days per week" | df_r_base_ctrl$att == "Always restricted", "4", ifelse(
        df_r_base_ctrl$att == "Baseline: No new requirements" | df_r_base_ctrl$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_ctrl$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_ctrl$att == "Baseline: No campaigns" | df_r_base_ctrl$att == "Occasional campaigns" | df_r_base_ctrl$att == "Frequent campaigns", "6", ifelse(
            df_r_base_ctrl$att == "Baseline: No increase of support" | df_r_base_ctrl$att == "Reducing prices by 15%" | df_r_base_ctrl$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_ctrl$policy = factor(df_r_base_ctrl$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))



################## Choice by country ###############
################## DE ###############
df_conj_DE <- subset(df_conj, country == "DE")
c_base_DE <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_DE, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_DE)


c_attribute1 <- c(0, c_base_DE$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_DE$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_DE$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_DE$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_DE$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_DE$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_DE$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_DE$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_DE$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_DE$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_DE$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_DE$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_DE$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_DE$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_DE$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_DE$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_DE$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_DE$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_DE$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_DE$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_DE$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_DE$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_DE$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_DE$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_DE$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_DE <- data.frame(c,s)
df_c_base_DE$att <- c(c_base_DE$attributes$attrib1lab[1], c_base_DE$attributes$attrib1lab[2:3], 
                      c_base_DE$attributes$attrib2lab[1], c_base_DE$attributes$attrib2lab[2:5],
                      c_base_DE$attributes$attrib7lab[3], c_base_DE$attributes$attrib7lab[2], c_base_DE$attributes$attrib7lab[1],
                      c_base_DE$attributes$attrib3lab[5], c_base_DE$attributes$attrib3lab[4], c_base_DE$attributes$attrib3lab[3], c_base_DE$attributes$attrib3lab[2], c_base_DE$attributes$attrib3lab[1], 
                      c_base_DE$attributes$attrib5lab[3], c_base_DE$attributes$attrib5lab[2], c_base_DE$attributes$attrib5lab[1],
                      c_base_DE$attributes$attrib4lab[3], c_base_DE$attributes$attrib4lab[2], c_base_DE$attributes$attrib4lab[1], 
                      c_base_DE$attributes$attrib6lab[3], c_base_DE$attributes$attrib6lab[2], c_base_DE$attributes$attrib6lab[1])
df_c_base_DE$att <- factor(df_c_base_DE$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_DE$policy <- ifelse(df_c_base_DE$att == "Baseline: No new tax" |  df_c_base_DE$att == "Increasing prices by 15%" |  df_c_base_DE$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_DE$att == "Baseline: General government budget" | df_c_base_DE$att == "Public environmental and climate protection programs" | df_c_base_DE$att == "Public programs for low-income households" | df_c_base_DE$att == "Reduce income taxes" | df_c_base_DE$att == "No tax revenues", "2", ifelse(
    df_c_base_DE$att == "Baseline: Keeping subsidies at current level" | df_c_base_DE$att == "Halving subsidies" | df_c_base_DE$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_DE$att == "Baseline: Never restricted" | df_c_base_DE$att == "1 day per week" | df_c_base_DE$att == "3 days per week" | df_c_base_DE$att == "5 days per week" | df_c_base_DE$att == "Always restricted", "4", ifelse(
        df_c_base_DE$att == "Baseline: No new requirements" | df_c_base_DE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_DE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_DE$att == "Baseline: No campaigns" | df_c_base_DE$att == "Occasional campaigns" | df_c_base_DE$att == "Frequent campaigns", "6", ifelse(
            df_c_base_DE$att == "Baseline: No increase of support" | df_c_base_DE$att == "Reducing prices by 15%" | df_c_base_DE$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_DE$policy = factor(df_c_base_DE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_c_base_DE <- ggplot(df_c_base_DE, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_c_base_DE

################## CN ###############

df_conj_CN <- subset(df_conj, country == "CN")
c_base_CN <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_CN, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_CN)


c_attribute1 <- c(0, c_base_CN$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_CN$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_CN$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_CN$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_CN$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_CN$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_CN$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_CN$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_CN$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_CN$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_CN$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_CN$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_CN$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_CN$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_CN$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_CN$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_CN$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_CN$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_CN$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_CN$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_CN$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_CN$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_CN$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_CN$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_CN$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_CN <- data.frame(c,s)
df_c_base_CN$att <- c(c_base_CN$attributes$attrib1lab[1], c_base_CN$attributes$attrib1lab[2:3], 
                      c_base_CN$attributes$attrib2lab[1], c_base_CN$attributes$attrib2lab[2:5],
                      c_base_CN$attributes$attrib7lab[3], c_base_CN$attributes$attrib7lab[2], c_base_CN$attributes$attrib7lab[1],
                      c_base_CN$attributes$attrib3lab[5], c_base_CN$attributes$attrib3lab[4], c_base_CN$attributes$attrib3lab[3], c_base_CN$attributes$attrib3lab[2], c_base_CN$attributes$attrib3lab[1], 
                      c_base_CN$attributes$attrib5lab[3], c_base_CN$attributes$attrib5lab[2], c_base_CN$attributes$attrib5lab[1],
                      c_base_CN$attributes$attrib4lab[3], c_base_CN$attributes$attrib4lab[2], c_base_CN$attributes$attrib4lab[1], 
                      c_base_CN$attributes$attrib6lab[3], c_base_CN$attributes$attrib6lab[2], c_base_CN$attributes$attrib6lab[1])
df_c_base_CN$att <- factor(df_c_base_CN$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_CN$policy <- ifelse(df_c_base_CN$att == "Baseline: No new tax" |  df_c_base_CN$att == "Increasing prices by 15%" |  df_c_base_CN$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_CN$att == "Baseline: General government budget" | df_c_base_CN$att == "Public environmental and climate protection programs" | df_c_base_CN$att == "Public programs for low-income households" | df_c_base_CN$att == "Reduce income taxes" | df_c_base_CN$att == "No tax revenues", "2", ifelse(
    df_c_base_CN$att == "Baseline: Keeping subsidies at current level" | df_c_base_CN$att == "Halving subsidies" | df_c_base_CN$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_CN$att == "Baseline: Never restricted" | df_c_base_CN$att == "1 day per week" | df_c_base_CN$att == "3 days per week" | df_c_base_CN$att == "5 days per week" | df_c_base_CN$att == "Always restricted", "4", ifelse(
        df_c_base_CN$att == "Baseline: No new requirements" | df_c_base_CN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_CN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_CN$att == "Baseline: No campaigns" | df_c_base_CN$att == "Occasional campaigns" | df_c_base_CN$att == "Frequent campaigns", "6", ifelse(
            df_c_base_CN$att == "Baseline: No increase of support" | df_c_base_CN$att == "Reducing prices by 15%" | df_c_base_CN$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_CN$policy = factor(df_c_base_CN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_c_base_CN <- ggplot(df_c_base_CN, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_c_base_CN

################## US ###############

df_conj_US <- subset(df_conj, country == "US")
c_base_US <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_US, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_US)


c_attribute1 <- c(0, c_base_US$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_US$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_US$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_US$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_US$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_US$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_US$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_US$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_US$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_US$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_US$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_US$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_US$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_US$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_US$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_US$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_US$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_US$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_US$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_US$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_US$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_US$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_US$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_US$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_US$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_US$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_US$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_US$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_US$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_US$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_US <- data.frame(c,s)
df_c_base_US$att <- c(c_base_US$attributes$attrib1lab[1], c_base_US$attributes$attrib1lab[2:3], 
                      c_base_US$attributes$attrib2lab[1], c_base_US$attributes$attrib2lab[2:5],
                      c_base_US$attributes$attrib7lab[3], c_base_US$attributes$attrib7lab[2], c_base_US$attributes$attrib7lab[1],
                      c_base_US$attributes$attrib3lab[5], c_base_US$attributes$attrib3lab[4], c_base_US$attributes$attrib3lab[3], c_base_US$attributes$attrib3lab[2], c_base_US$attributes$attrib3lab[1], 
                      c_base_US$attributes$attrib5lab[3], c_base_US$attributes$attrib5lab[2], c_base_US$attributes$attrib5lab[1],
                      c_base_US$attributes$attrib4lab[3], c_base_US$attributes$attrib4lab[2], c_base_US$attributes$attrib4lab[1], 
                      c_base_US$attributes$attrib6lab[3], c_base_US$attributes$attrib6lab[2], c_base_US$attributes$attrib6lab[1])
df_c_base_US$att <- factor(df_c_base_US$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_US$policy <- ifelse(df_c_base_US$att == "Baseline: No new tax" |  df_c_base_US$att == "Increasing prices by 15%" |  df_c_base_US$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_US$att == "Baseline: General government budget" | df_c_base_US$att == "Public environmental and climate protection programs" | df_c_base_US$att == "Public programs for low-income households" | df_c_base_US$att == "Reduce income taxes" | df_c_base_US$att == "No tax revenues", "2", ifelse(
    df_c_base_US$att == "Baseline: Keeping subsidies at current level" | df_c_base_US$att == "Halving subsidies" | df_c_base_US$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_US$att == "Baseline: Never restricted" | df_c_base_US$att == "1 day per week" | df_c_base_US$att == "3 days per week" | df_c_base_US$att == "5 days per week" | df_c_base_US$att == "Always restricted", "4", ifelse(
        df_c_base_US$att == "Baseline: No new requirements" | df_c_base_US$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_US$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_US$att == "Baseline: No campaigns" | df_c_base_US$att == "Occasional campaigns" | df_c_base_US$att == "Frequent campaigns", "6", ifelse(
            df_c_base_US$att == "Baseline: No increase of support" | df_c_base_US$att == "Reducing prices by 15%" | df_c_base_US$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_US$policy = factor(df_c_base_US$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_c_base_US <- ggplot(df_c_base_US, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_c_base_US

################## Rating by country ###############
################## DE ###############
df_conj_DE <- subset(df_conj, country == "DE")
r_base_DE <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_DE, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_DE)


c_attribute1 <- c(0, r_base_DE$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_DE$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_DE$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_DE$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_DE$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_DE$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_DE$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_DE$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_DE$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_DE$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_DE$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_DE$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_DE$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_DE$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_DE$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_DE$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_DE$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_DE$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_DE$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_DE$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_DE$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_DE$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_DE$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_DE$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_DE$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_DE <- data.frame(c,s)
df_r_base_DE$att <- c(r_base_DE$attributes$attrib1lab[1], r_base_DE$attributes$attrib1lab[2:3], 
                      r_base_DE$attributes$attrib2lab[1], r_base_DE$attributes$attrib2lab[2:5],
                      r_base_DE$attributes$attrib7lab[3], r_base_DE$attributes$attrib7lab[2], r_base_DE$attributes$attrib7lab[1],
                      r_base_DE$attributes$attrib3lab[5], r_base_DE$attributes$attrib3lab[4], r_base_DE$attributes$attrib3lab[3], r_base_DE$attributes$attrib3lab[2], r_base_DE$attributes$attrib3lab[1], 
                      r_base_DE$attributes$attrib5lab[3], r_base_DE$attributes$attrib5lab[2], r_base_DE$attributes$attrib5lab[1],
                      r_base_DE$attributes$attrib4lab[3], r_base_DE$attributes$attrib4lab[2], r_base_DE$attributes$attrib4lab[1], 
                      r_base_DE$attributes$attrib6lab[3], r_base_DE$attributes$attrib6lab[2], r_base_DE$attributes$attrib6lab[1])
df_r_base_DE$att <- factor(df_r_base_DE$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_DE$policy <- ifelse(df_r_base_DE$att == "Baseline: No new tax" |  df_r_base_DE$att == "Increasing prices by 15%" |  df_r_base_DE$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_DE$att == "Baseline: General government budget" | df_r_base_DE$att == "Public environmental and climate protection programs" | df_r_base_DE$att == "Public programs for low-income households" | df_r_base_DE$att == "Reduce income taxes" | df_r_base_DE$att == "No tax revenues", "2", ifelse(
    df_r_base_DE$att == "Baseline: Keeping subsidies at current level" | df_r_base_DE$att == "Halving subsidies" | df_r_base_DE$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_DE$att == "Baseline: Never restricted" | df_r_base_DE$att == "1 day per week" | df_r_base_DE$att == "3 days per week" | df_r_base_DE$att == "5 days per week" | df_r_base_DE$att == "Always restricted", "4", ifelse(
        df_r_base_DE$att == "Baseline: No new requirements" | df_r_base_DE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_DE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_DE$att == "Baseline: No campaigns" | df_r_base_DE$att == "Occasional campaigns" | df_r_base_DE$att == "Frequent campaigns", "6", ifelse(
            df_r_base_DE$att == "Baseline: No increase of support" | df_r_base_DE$att == "Reducing prices by 15%" | df_r_base_DE$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_DE$policy = factor(df_r_base_DE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_r_base_DE <- ggplot(df_r_base_DE, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_r_base_DE

################## CN ###############

df_conj_CN <- subset(df_conj, country == "CN")
r_base_CN <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_CN, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_CN)


c_attribute1 <- c(0, r_base_CN$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_CN$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_CN$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_CN$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_CN$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_CN$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_CN$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_CN$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_CN$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_CN$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_CN$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_CN$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_CN$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_CN$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_CN$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_CN$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_CN$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_CN$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_CN$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_CN$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_CN$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_CN$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_CN$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_CN$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_CN$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_CN <- data.frame(c,s)
df_r_base_CN$att <- c(r_base_CN$attributes$attrib1lab[1], r_base_CN$attributes$attrib1lab[2:3], 
                      r_base_CN$attributes$attrib2lab[1], r_base_CN$attributes$attrib2lab[2:5],
                      r_base_CN$attributes$attrib7lab[3], r_base_CN$attributes$attrib7lab[2], r_base_CN$attributes$attrib7lab[1],
                      r_base_CN$attributes$attrib3lab[5], r_base_CN$attributes$attrib3lab[4], r_base_CN$attributes$attrib3lab[3], r_base_CN$attributes$attrib3lab[2], r_base_CN$attributes$attrib3lab[1], 
                      r_base_CN$attributes$attrib5lab[3], r_base_CN$attributes$attrib5lab[2], r_base_CN$attributes$attrib5lab[1],
                      r_base_CN$attributes$attrib4lab[3], r_base_CN$attributes$attrib4lab[2], r_base_CN$attributes$attrib4lab[1], 
                      r_base_CN$attributes$attrib6lab[3], r_base_CN$attributes$attrib6lab[2], r_base_CN$attributes$attrib6lab[1])
df_r_base_CN$att <- factor(df_r_base_CN$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_CN$policy <- ifelse(df_r_base_CN$att == "Baseline: No new tax" |  df_r_base_CN$att == "Increasing prices by 15%" |  df_r_base_CN$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_CN$att == "Baseline: General government budget" | df_r_base_CN$att == "Public environmental and climate protection programs" | df_r_base_CN$att == "Public programs for low-income households" | df_r_base_CN$att == "Reduce income taxes" | df_r_base_CN$att == "No tax revenues", "2", ifelse(
    df_r_base_CN$att == "Baseline: Keeping subsidies at current level" | df_r_base_CN$att == "Halving subsidies" | df_r_base_CN$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_CN$att == "Baseline: Never restricted" | df_r_base_CN$att == "1 day per week" | df_r_base_CN$att == "3 days per week" | df_r_base_CN$att == "5 days per week" | df_r_base_CN$att == "Always restricted", "4", ifelse(
        df_r_base_CN$att == "Baseline: No new requirements" | df_r_base_CN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_CN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_CN$att == "Baseline: No campaigns" | df_r_base_CN$att == "Occasional campaigns" | df_r_base_CN$att == "Frequent campaigns", "6", ifelse(
            df_r_base_CN$att == "Baseline: No increase of support" | df_r_base_CN$att == "Reducing prices by 15%" | df_r_base_CN$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_CN$policy = factor(df_r_base_CN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_r_base_CN <- ggplot(df_r_base_CN, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_r_base_CN


################## US ###############

df_conj_US <- subset(df_conj, country == "US")
r_base_US <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_US, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_US)


c_attribute1 <- c(0, r_base_US$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_US$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_US$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_US$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_US$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_US$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_US$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_US$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_US$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_US$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_US$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_US$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_US$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_US$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_US$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_US$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_US$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_US$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_US$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_US$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_US$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_US$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_US$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_US$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_US$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_US$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_US$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_US$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_US$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_US$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_US <- data.frame(c,s)
df_r_base_US$att <- c(r_base_US$attributes$attrib1lab[1], r_base_US$attributes$attrib1lab[2:3], 
                      r_base_US$attributes$attrib2lab[1], r_base_US$attributes$attrib2lab[2:5],
                      r_base_US$attributes$attrib7lab[3], r_base_US$attributes$attrib7lab[2], r_base_US$attributes$attrib7lab[1],
                      r_base_US$attributes$attrib3lab[5], r_base_US$attributes$attrib3lab[4], r_base_US$attributes$attrib3lab[3], r_base_US$attributes$attrib3lab[2], r_base_US$attributes$attrib3lab[1], 
                      r_base_US$attributes$attrib5lab[3], r_base_US$attributes$attrib5lab[2], r_base_US$attributes$attrib5lab[1],
                      r_base_US$attributes$attrib4lab[3], r_base_US$attributes$attrib4lab[2], r_base_US$attributes$attrib4lab[1], 
                      r_base_US$attributes$attrib6lab[3], r_base_US$attributes$attrib6lab[2], r_base_US$attributes$attrib6lab[1])
df_r_base_US$att <- factor(df_r_base_US$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_US$policy <- ifelse(df_r_base_US$att == "Baseline: No new tax" |  df_r_base_US$att == "Increasing prices by 15%" |  df_r_base_US$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_US$att == "Baseline: General government budget" | df_r_base_US$att == "Public environmental and climate protection programs" | df_r_base_US$att == "Public programs for low-income households" | df_r_base_US$att == "Reduce income taxes" | df_r_base_US$att == "No tax revenues", "2", ifelse(
    df_r_base_US$att == "Baseline: Keeping subsidies at current level" | df_r_base_US$att == "Halving subsidies" | df_r_base_US$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_US$att == "Baseline: Never restricted" | df_r_base_US$att == "1 day per week" | df_r_base_US$att == "3 days per week" | df_r_base_US$att == "5 days per week" | df_r_base_US$att == "Always restricted", "4", ifelse(
        df_r_base_US$att == "Baseline: No new requirements" | df_r_base_US$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_US$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_US$att == "Baseline: No campaigns" | df_r_base_US$att == "Occasional campaigns" | df_r_base_US$att == "Frequent campaigns", "6", ifelse(
            df_r_base_US$att == "Baseline: No increase of support" | df_r_base_US$att == "Reducing prices by 15%" | df_r_base_US$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_US$policy = factor(df_r_base_US$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_r_base_US <- ggplot(df_r_base_US, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  facet_grid(policy ~ ., scales = "free")
p_r_base_US

################## Marginal Means ###############

library(cregg)


#choice
cmm <- mm(df_conj, choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
            attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(cmm)

cmm$att <- factor(cmm$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                              "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                              "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                              "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                              "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                              "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                              "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                  labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                   "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                   "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                   "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                   "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                   "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                   "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

cmm$poltyp <- factor(ifelse(cmm$att == "No new tax" |  cmm$att == "Increasing prices by 15%" |  cmm$att == "Increasing prices by 30%" | 
                              cmm$att == "General government budget" | cmm$att == "Public environmental and climate protection programs" | cmm$att == "Public programs for low-income households" | cmm$att == "Reduce income taxes" | cmm$att == "No tax revenues" |
                              cmm$att == "Keeping subsidies at current level" | cmm$att == "Halving subsidies" | cmm$att == "Eliminating subsidies",
                            "Market-based / push", ifelse(
                              cmm$att == "Never restricted" | cmm$att == "1 day per week" | cmm$att == "3 days per week" | cmm$att == "5 days per week" | cmm$att == "Always restricted" | 
                                cmm$att == "No new requirements" | cmm$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmm$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                              "Non-market-based / push", ifelse(
                                cmm$att == "No campaigns" | cmm$att == "Occasional campaigns" | cmm$att == "Frequent campaigns", 
                                "Non-market-based / pull", ifelse(
                                  cmm$att == "No increase of support" | cmm$att == "Reducing prices by 15%" | cmm$att == "Reducing prices by 30%",
                                  "Market-based / pull", "Other Policy Features")))))

cmm$policy <- ifelse(cmm$att == "No new tax" |  cmm$att == "Increasing prices by 15%" |  cmm$att == "Increasing prices by 30%", "1", ifelse(
  cmm$att == "General government budget" | cmm$att == "Public environmental and climate protection programs" | cmm$att == "Public programs for low-income households" | cmm$att == "Reduce income taxes" | cmm$att == "No tax revenues", "2", ifelse(
    cmm$att == "Keeping subsidies at current level" | cmm$att == "Halving subsidies" | cmm$att == "Eliminating subsidies", "3", ifelse(
      cmm$att == "Never restricted" | cmm$att == "1 day per week" | cmm$att == "3 days per week" | cmm$att == "5 days per week" | cmm$att == "Always restricted", "4", ifelse(
        cmm$att == "No new requirements" | cmm$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmm$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          cmm$att == "No campaigns" | cmm$att == "Occasional campaigns" | cmm$att == "Frequent campaigns", "6", ifelse(
            cmm$att == "No increase of support" | cmm$att == "Reducing prices by 15%" | cmm$att == "Reducing prices by 30%", "7",
            "8")))))))

cmm$policy = factor(cmm$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_cmm <- ggplot(cmm, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_cmm


#rate
rmm <- mm(df_conj, rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
            attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(rmm)

rmm$att <- factor(rmm$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                              "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                              "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                              "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                              "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                              "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                              "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                  labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                   "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                   "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                   "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                   "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                   "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                   "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

rmm$poltyp <- factor(ifelse(rmm$att == "No new tax" |  rmm$att == "Increasing prices by 15%" |  rmm$att == "Increasing prices by 30%" | 
                              rmm$att == "General government budget" | rmm$att == "Public environmental and climate protection programs" | rmm$att == "Public programs for low-income households" | rmm$att == "Reduce income taxes" | rmm$att == "No tax revenues" |
                              rmm$att == "Keeping subsidies at current level" | rmm$att == "Halving subsidies" | rmm$att == "Eliminating subsidies",
                            "Market-based / push", ifelse(
                              rmm$att == "Never restricted" | rmm$att == "1 day per week" | rmm$att == "3 days per week" | rmm$att == "5 days per week" | rmm$att == "Always restricted" | 
                                rmm$att == "No new requirements" | rmm$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmm$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                              "Non-market-based / push", ifelse(
                                rmm$att == "No campaigns" | rmm$att == "Occasional campaigns" | rmm$att == "Frequent campaigns", 
                                "Non-market-based / pull", ifelse(
                                  rmm$att == "No increase of support" | rmm$att == "Reducing prices by 15%" | rmm$att == "Reducing prices by 30%",
                                  "Market-based / pull", "Other Policy Features")))))

rmm$policy <- ifelse(rmm$att == "No new tax" |  rmm$att == "Increasing prices by 15%" |  rmm$att == "Increasing prices by 30%", "1", ifelse(
  rmm$att == "General government budget" | rmm$att == "Public environmental and climate protection programs" | rmm$att == "Public programs for low-income households" | rmm$att == "Reduce income taxes" | rmm$att == "No tax revenues", "2", ifelse(
    rmm$att == "Keeping subsidies at current level" | rmm$att == "Halving subsidies" | rmm$att == "Eliminating subsidies", "3", ifelse(
      rmm$att == "Never restricted" | rmm$att == "1 day per week" | rmm$att == "3 days per week" | rmm$att == "5 days per week" | rmm$att == "Always restricted", "4", ifelse(
        rmm$att == "No new requirements" | rmm$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmm$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          rmm$att == "No campaigns" | rmm$att == "Occasional campaigns" | rmm$att == "Frequent campaigns", "6", ifelse(
            rmm$att == "No increase of support" | rmm$att == "Reducing prices by 15%" | rmm$att == "Reducing prices by 30%", "7",
            "8")))))))

rmm$policy = factor(rmm$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_rmm <- ggplot(rmm, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_rmm



################## Marginal Means by country ###############
#choice DE
cmmDE <- mm(df_conj_DE, choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(cmmDE)

cmmDE$att <- factor(cmmDE$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

cmmDE$poltyp <- factor(ifelse(cmmDE$att == "No new tax" |  cmmDE$att == "Increasing prices by 15%" |  cmmDE$att == "Increasing prices by 30%" | 
                                cmmDE$att == "General government budget" | cmmDE$att == "Public environmental and climate protection programs" | cmmDE$att == "Public programs for low-income households" | cmmDE$att == "Reduce income taxes" | cmmDE$att == "No tax revenues" |
                                cmmDE$att == "Keeping subsidies at current level" | cmmDE$att == "Halving subsidies" | cmmDE$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                cmmDE$att == "Never restricted" | cmmDE$att == "1 day per week" | cmmDE$att == "3 days per week" | cmmDE$att == "5 days per week" | cmmDE$att == "Always restricted" | 
                                  cmmDE$att == "No new requirements" | cmmDE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmDE$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  cmmDE$att == "No campaigns" | cmmDE$att == "Occasional campaigns" | cmmDE$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    cmmDE$att == "No increase of support" | cmmDE$att == "Reducing prices by 15%" | cmmDE$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

cmmDE$policy <- ifelse(cmmDE$att == "No new tax" |  cmmDE$att == "Increasing prices by 15%" |  cmmDE$att == "Increasing prices by 30%", "1", ifelse(
  cmmDE$att == "General government budget" | cmmDE$att == "Public environmental and climate protection programs" | cmmDE$att == "Public programs for low-income households" | cmmDE$att == "Reduce income taxes" | cmmDE$att == "No tax revenues", "2", ifelse(
    cmmDE$att == "Keeping subsidies at current level" | cmmDE$att == "Halving subsidies" | cmmDE$att == "Eliminating subsidies", "3", ifelse(
      cmmDE$att == "Never restricted" | cmmDE$att == "1 day per week" | cmmDE$att == "3 days per week" | cmmDE$att == "5 days per week" | cmmDE$att == "Always restricted", "4", ifelse(
        cmmDE$att == "No new requirements" | cmmDE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmDE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          cmmDE$att == "No campaigns" | cmmDE$att == "Occasional campaigns" | cmmDE$att == "Frequent campaigns", "6", ifelse(
            cmmDE$att == "No increase of support" | cmmDE$att == "Reducing prices by 15%" | cmmDE$att == "Reducing prices by 30%", "7",
            "8")))))))

cmmDE$policy = factor(cmmDE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_cmmDE <- ggplot(cmmDE, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_cmmDE


#rate DE
rmmDE <- mm(df_conj_DE, rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(rmmDE)

rmmDE$att <- factor(rmmDE$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

rmmDE$poltyp <- factor(ifelse(rmmDE$att == "No new tax" |  rmmDE$att == "Increasing prices by 15%" |  rmmDE$att == "Increasing prices by 30%" | 
                                rmmDE$att == "General government budget" | rmmDE$att == "Public environmental and climate protection programs" | rmmDE$att == "Public programs for low-income households" | rmmDE$att == "Reduce income taxes" | rmmDE$att == "No tax revenues" |
                                rmmDE$att == "Keeping subsidies at current level" | rmmDE$att == "Halving subsidies" | rmmDE$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                rmmDE$att == "Never restricted" | rmmDE$att == "1 day per week" | rmmDE$att == "3 days per week" | rmmDE$att == "5 days per week" | rmmDE$att == "Always restricted" | 
                                  rmmDE$att == "No new requirements" | rmmDE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmDE$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  rmmDE$att == "No campaigns" | rmmDE$att == "Occasional campaigns" | rmmDE$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    rmmDE$att == "No increase of support" | rmmDE$att == "Reducing prices by 15%" | rmmDE$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

rmmDE$policy <- ifelse(rmmDE$att == "No new tax" |  rmmDE$att == "Increasing prices by 15%" |  rmmDE$att == "Increasing prices by 30%", "1", ifelse(
  rmmDE$att == "General government budget" | rmmDE$att == "Public environmental and climate protection programs" | rmmDE$att == "Public programs for low-income households" | rmmDE$att == "Reduce income taxes" | rmmDE$att == "No tax revenues", "2", ifelse(
    rmmDE$att == "Keeping subsidies at current level" | rmmDE$att == "Halving subsidies" | rmmDE$att == "Eliminating subsidies", "3", ifelse(
      rmmDE$att == "Never restricted" | rmmDE$att == "1 day per week" | rmmDE$att == "3 days per week" | rmmDE$att == "5 days per week" | rmmDE$att == "Always restricted", "4", ifelse(
        rmmDE$att == "No new requirements" | rmmDE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmDE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          rmmDE$att == "No campaigns" | rmmDE$att == "Occasional campaigns" | rmmDE$att == "Frequent campaigns", "6", ifelse(
            rmmDE$att == "No increase of support" | rmmDE$att == "Reducing prices by 15%" | rmmDE$att == "Reducing prices by 30%", "7",
            "8")))))))

rmmDE$policy = factor(rmmDE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_rmmDE <- ggplot(rmmDE, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_rmmDE

#choice CN
cmmCN <- mm(df_conj_CN, choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(cmmCN)

cmmCN$att <- factor(cmmCN$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

cmmCN$poltyp <- factor(ifelse(cmmCN$att == "No new tax" |  cmmCN$att == "Increasing prices by 15%" |  cmmCN$att == "Increasing prices by 30%" | 
                                cmmCN$att == "General government budget" | cmmCN$att == "Public environmental and climate protection programs" | cmmCN$att == "Public programs for low-income households" | cmmCN$att == "Reduce income taxes" | cmmCN$att == "No tax revenues" |
                                cmmCN$att == "Keeping subsidies at current level" | cmmCN$att == "Halving subsidies" | cmmCN$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                cmmCN$att == "Never restricted" | cmmCN$att == "1 day per week" | cmmCN$att == "3 days per week" | cmmCN$att == "5 days per week" | cmmCN$att == "Always restricted" | 
                                  cmmCN$att == "No new requirements" | cmmCN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmCN$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  cmmCN$att == "No campaigns" | cmmCN$att == "Occasional campaigns" | cmmCN$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    cmmCN$att == "No increase of support" | cmmCN$att == "Reducing prices by 15%" | cmmCN$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

cmmCN$policy <- ifelse(cmmCN$att == "No new tax" |  cmmCN$att == "Increasing prices by 15%" |  cmmCN$att == "Increasing prices by 30%", "1", ifelse(
  cmmCN$att == "General government budget" | cmmCN$att == "Public environmental and climate protection programs" | cmmCN$att == "Public programs for low-income households" | cmmCN$att == "Reduce income taxes" | cmmCN$att == "No tax revenues", "2", ifelse(
    cmmCN$att == "Keeping subsidies at current level" | cmmCN$att == "Halving subsidies" | cmmCN$att == "Eliminating subsidies", "3", ifelse(
      cmmCN$att == "Never restricted" | cmmCN$att == "1 day per week" | cmmCN$att == "3 days per week" | cmmCN$att == "5 days per week" | cmmCN$att == "Always restricted", "4", ifelse(
        cmmCN$att == "No new requirements" | cmmCN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmCN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          cmmCN$att == "No campaigns" | cmmCN$att == "Occasional campaigns" | cmmCN$att == "Frequent campaigns", "6", ifelse(
            cmmCN$att == "No increase of support" | cmmCN$att == "Reducing prices by 15%" | cmmCN$att == "Reducing prices by 30%", "7",
            "8")))))))

cmmCN$policy = factor(cmmCN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_cmmCN <- ggplot(cmmCN, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_cmmCN


#rate CN
rmmCN <- mm(df_conj_CN, rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(rmmCN)

rmmCN$att <- factor(rmmCN$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

rmmCN$poltyp <- factor(ifelse(rmmCN$att == "No new tax" |  rmmCN$att == "Increasing prices by 15%" |  rmmCN$att == "Increasing prices by 30%" | 
                                rmmCN$att == "General government budget" | rmmCN$att == "Public environmental and climate protection programs" | rmmCN$att == "Public programs for low-income households" | rmmCN$att == "Reduce income taxes" | rmmCN$att == "No tax revenues" |
                                rmmCN$att == "Keeping subsidies at current level" | rmmCN$att == "Halving subsidies" | rmmCN$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                rmmCN$att == "Never restricted" | rmmCN$att == "1 day per week" | rmmCN$att == "3 days per week" | rmmCN$att == "5 days per week" | rmmCN$att == "Always restricted" | 
                                  rmmCN$att == "No new requirements" | rmmCN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmCN$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  rmmCN$att == "No campaigns" | rmmCN$att == "Occasional campaigns" | rmmCN$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    rmmCN$att == "No increase of support" | rmmCN$att == "Reducing prices by 15%" | rmmCN$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

rmmCN$policy <- ifelse(rmmCN$att == "No new tax" |  rmmCN$att == "Increasing prices by 15%" |  rmmCN$att == "Increasing prices by 30%", "1", ifelse(
  rmmCN$att == "General government budget" | rmmCN$att == "Public environmental and climate protection programs" | rmmCN$att == "Public programs for low-income households" | rmmCN$att == "Reduce income taxes" | rmmCN$att == "No tax revenues", "2", ifelse(
    rmmCN$att == "Keeping subsidies at current level" | rmmCN$att == "Halving subsidies" | rmmCN$att == "Eliminating subsidies", "3", ifelse(
      rmmCN$att == "Never restricted" | rmmCN$att == "1 day per week" | rmmCN$att == "3 days per week" | rmmCN$att == "5 days per week" | rmmCN$att == "Always restricted", "4", ifelse(
        rmmCN$att == "No new requirements" | rmmCN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmCN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          rmmCN$att == "No campaigns" | rmmCN$att == "Occasional campaigns" | rmmCN$att == "Frequent campaigns", "6", ifelse(
            rmmCN$att == "No increase of support" | rmmCN$att == "Reducing prices by 15%" | rmmCN$att == "Reducing prices by 30%", "7",
            "8")))))))

rmmCN$policy = factor(rmmCN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_rmmCN <- ggplot(rmmCN, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_rmmCN


#choice US
cmmUS <- mm(df_conj_US, choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(cmmUS)

cmmUS$att <- factor(cmmUS$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

cmmUS$poltyp <- factor(ifelse(cmmUS$att == "No new tax" |  cmmUS$att == "Increasing prices by 15%" |  cmmUS$att == "Increasing prices by 30%" | 
                                cmmUS$att == "General government budget" | cmmUS$att == "Public environmental and climate protection programs" | cmmUS$att == "Public programs for low-income households" | cmmUS$att == "Reduce income taxes" | cmmUS$att == "No tax revenues" |
                                cmmUS$att == "Keeping subsidies at current level" | cmmUS$att == "Halving subsidies" | cmmUS$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                cmmUS$att == "Never restricted" | cmmUS$att == "1 day per week" | cmmUS$att == "3 days per week" | cmmUS$att == "5 days per week" | cmmUS$att == "Always restricted" | 
                                  cmmUS$att == "No new requirements" | cmmUS$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmUS$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  cmmUS$att == "No campaigns" | cmmUS$att == "Occasional campaigns" | cmmUS$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    cmmUS$att == "No increase of support" | cmmUS$att == "Reducing prices by 15%" | cmmUS$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

cmmUS$policy <- ifelse(cmmUS$att == "No new tax" |  cmmUS$att == "Increasing prices by 15%" |  cmmUS$att == "Increasing prices by 30%", "1", ifelse(
  cmmUS$att == "General government budget" | cmmUS$att == "Public environmental and climate protection programs" | cmmUS$att == "Public programs for low-income households" | cmmUS$att == "Reduce income taxes" | cmmUS$att == "No tax revenues", "2", ifelse(
    cmmUS$att == "Keeping subsidies at current level" | cmmUS$att == "Halving subsidies" | cmmUS$att == "Eliminating subsidies", "3", ifelse(
      cmmUS$att == "Never restricted" | cmmUS$att == "1 day per week" | cmmUS$att == "3 days per week" | cmmUS$att == "5 days per week" | cmmUS$att == "Always restricted", "4", ifelse(
        cmmUS$att == "No new requirements" | cmmUS$att == "Maximum consumption of 5.9 litres per 100 kilometres" | cmmUS$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          cmmUS$att == "No campaigns" | cmmUS$att == "Occasional campaigns" | cmmUS$att == "Frequent campaigns", "6", ifelse(
            cmmUS$att == "No increase of support" | cmmUS$att == "Reducing prices by 15%" | cmmUS$att == "Reducing prices by 30%", "7",
            "8")))))))

cmmUS$policy = factor(cmmUS$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_cmmUS <- ggplot(cmmUS, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_cmmUS


#rate US
rmmUS <- mm(df_conj_US, rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab + attrib5_lab +  
              attrib4_lab + attrib6_lab, id = ~id, design=design, baselines = baselines_base)

plot(rmmUS)

rmmUS$att <- factor(rmmUS$level, levels = c(rev(c("3) No new tax", "2) Increasing prices by 15%", "1) Increasing prices by 30%",
                                                  "2) General government budget", "1) Public environmental and climate protection programs", "3) Public programs for low-income households", "4) Reduce income taxes", "5) No tax revenues",
                                                  "3) Keeping subsidies at current level", "2) Halving subsidies ", "1) Eliminating subsidies",
                                                  "5) Never restricted", "4) 1 day per week", "3) 3 days per week", "2) 5 days per week", "1) Always restricted", 
                                                  "3) No new requirements", "2) At least 40 miles per gallon ", "1) At least 55 miles per gallon ",
                                                  "3) No increase of support", "2) Reducing prices by 15%", "1) Reducing prices by 30%",
                                                  "3) No campaigns", "2) Occasional campaigns", "1) Frequent campaigns"))),
                    labels = c(rev(c("No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                     "General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                     "Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                     "Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                     "No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                     "No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                     "No campaigns", "Occasional campaigns", "Frequent campaigns"))))

rmmUS$poltyp <- factor(ifelse(rmmUS$att == "No new tax" |  rmmUS$att == "Increasing prices by 15%" |  rmmUS$att == "Increasing prices by 30%" | 
                                rmmUS$att == "General government budget" | rmmUS$att == "Public environmental and climate protection programs" | rmmUS$att == "Public programs for low-income households" | rmmUS$att == "Reduce income taxes" | rmmUS$att == "No tax revenues" |
                                rmmUS$att == "Keeping subsidies at current level" | rmmUS$att == "Halving subsidies" | rmmUS$att == "Eliminating subsidies",
                              "Market-based / push", ifelse(
                                rmmUS$att == "Never restricted" | rmmUS$att == "1 day per week" | rmmUS$att == "3 days per week" | rmmUS$att == "5 days per week" | rmmUS$att == "Always restricted" | 
                                  rmmUS$att == "No new requirements" | rmmUS$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmUS$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                "Non-market-based / push", ifelse(
                                  rmmUS$att == "No campaigns" | rmmUS$att == "Occasional campaigns" | rmmUS$att == "Frequent campaigns", 
                                  "Non-market-based / pull", ifelse(
                                    rmmUS$att == "No increase of support" | rmmUS$att == "Reducing prices by 15%" | rmmUS$att == "Reducing prices by 30%",
                                    "Market-based / pull", "Other Policy Features")))))

rmmUS$policy <- ifelse(rmmUS$att == "No new tax" |  rmmUS$att == "Increasing prices by 15%" |  rmmUS$att == "Increasing prices by 30%", "1", ifelse(
  rmmUS$att == "General government budget" | rmmUS$att == "Public environmental and climate protection programs" | rmmUS$att == "Public programs for low-income households" | rmmUS$att == "Reduce income taxes" | rmmUS$att == "No tax revenues", "2", ifelse(
    rmmUS$att == "Keeping subsidies at current level" | rmmUS$att == "Halving subsidies" | rmmUS$att == "Eliminating subsidies", "3", ifelse(
      rmmUS$att == "Never restricted" | rmmUS$att == "1 day per week" | rmmUS$att == "3 days per week" | rmmUS$att == "5 days per week" | rmmUS$att == "Always restricted", "4", ifelse(
        rmmUS$att == "No new requirements" | rmmUS$att == "Maximum consumption of 5.9 litres per 100 kilometres" | rmmUS$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          rmmUS$att == "No campaigns" | rmmUS$att == "Occasional campaigns" | rmmUS$att == "Frequent campaigns", "6", ifelse(
            rmmUS$att == "No increase of support" | rmmUS$att == "Reducing prices by 15%" | rmmUS$att == "Reducing prices by 30%", "7",
            "8")))))))

rmmUS$policy = factor(rmmUS$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_rmmUS <- ggplot(rmmUS, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  facet_grid(policy ~ ., scales = "free")
p_rmmUS

############################# Rating overview (with 95% confidence interval) ###############################
library(Rmisc)
library(ggplot2)
library(cowplot)
library(dplyr)

#library(effects)
#library(lme4)
#library(data.table)

#library(nFactors)
#library(texreg)
#library(tidyr)
#library(forcats)
#library(scales)
#library(Matrix)

#country
summary(df_conj$rate[df_conj$country=="CN"])
summary(df_conj$rate[df_conj$country=="DE"])
summary(df_conj$rate[df_conj$country=="US"])
rate_summary_country <- summarySE(df_conj, measurevar="rate", groupvars=c("country"), na.rm = T, conf.interval = 0.95)
rate_summary_country

fig_rate_summary_country1 <- ggplot(rate_summary_country, aes(x=country, y=rate, fill=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=0.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab(" ") +
  ylab("Average rating") +
  #theme(plot.title = element_text(size = 20, face = "bold"),axis.text.x =element_text(size=20, face ="bold"),axis.text.y =element_text(size=20, face ="bold"),strip.text.x = element_text(size=20, face ="bold"), axis.title.y = element_text(size = 20, face ="bold"))+
  #ggtitle("Average Rating (1-7) by Country") +
  scale_fill_manual("", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("China", "Germany", "USA")) +
  coord_cartesian(ylim = c(1.27, 6.73)) 

legend3 <- get_legend(fig_rate_summary_country1)

fig_rate_summary_country <- fig_rate_summary_country1 + guides(fill=FALSE)

#overview average rating
rateavrg <- aggregate(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                        attrib5_lab +  attrib6_lab +  attrib7_lab + country, df_conj, mean)

fig_rate_average <- ggplot(rateavrg, aes(x=rate, colour=country)) + 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Distribution of average ratings for policy-packages") +
  ylab("Density") +
  scale_colour_manual(values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  coord_cartesian(ylim = c(.0257, 0.55), xlim = c(1.27, 6.73)) +
  #geom_vline(xintercept = 4.740764, color = "#e66101", size=1.2, linetype = "longdash") +
  #geom_vline(xintercept = 3.651830, color = "#5e3c99", size=1.2, linetype = "longdash") +
  #geom_vline(xintercept = 4.069573, color = "#b2abd2", size=1.2, linetype = "longdash") +
  guides(colour=FALSE) +
  geom_vline(xintercept = 4, lty="dashed")

#age
rate_summary_age <- summarySE(df_conj, measurevar="rate", groupvars=c("age_grp", "country"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_age2 <- ggplot(rate_summary_age, aes(x=age_grp, y=rate, fill=country, shape=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.5,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("") +
  ylab("Rating") +
  #theme(plot.title = element_text(size = 20, face = "bold"),axis.text.x =element_text(size=20, face ="bold"),axis.text.y =element_text(size=20, face ="bold"),strip.text.x = element_text(size=20, face ="bold"), axis.title.y = element_text(size = 20, face ="bold"))+
  ggtitle("Average rating by age groups") +
  scale_fill_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_shape_discrete() +
  coord_cartesian(ylim = c(1.27, 6.73)) +
  guides(fill=FALSE)


#gender
rate_summary_gender <- summarySE(df_conj, measurevar="rate", groupvars=c("gender", "country"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_gender1 <- ggplot(rate_summary_gender, aes(x=gender, y=rate, fill=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("") +
  ylab("Rating") +
  #theme(plot.title = element_text(size = 14, face = "bold"),axis.text.x =element_text(size=14, face ="bold"),axis.text.y =element_text(size=14, face ="bold"),strip.text.x = element_text(size=14, face ="bold"), axis.title.y = element_text(size = 14, face ="bold"))+
  ggtitle("Average rating by gender") +
  scale_fill_manual(name="Country", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("Female", "Male")) +
  coord_cartesian(ylim = c(1.27, 6.73))

legend2 <- get_legend(fig_rate_summary_gender1)

fig_rate_summary_gender <- fig_rate_summary_gender1 + guides(fill=FALSE)


#income
rate_summary_incgrp <- summarySE(df_conj, measurevar="rate", groupvars=c("incgrp", "country"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_incgrp <- ggplot(data=subset(rate_summary_incgrp, !is.na(incgrp)), aes(x=incgrp, y=rate, fill=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("") +
  ylab("Rating") +
  #theme(plot.title = element_text(size = 14, face = "bold"),axis.text.x =element_text(size=14, face ="bold"),axis.text.y =element_text(size=14, face ="bold"),strip.text.x = element_text(size=14, face ="bold"), axis.title.y = element_text(size = 14, face ="bold"))+
  ggtitle("Average rating by income groups (quintiles)") +
  scale_fill_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("1st Quintile", "2nd Quintile", "3rd Quintile", "4th Quintile", "5th Quintile")) +
  coord_cartesian(ylim = c(1.27, 6.73)) +
  guides(fill=FALSE)
fig_rate_summary_incgrp


#educ
rate_summary_educ <- summarySE(df_conj, measurevar="rate", groupvars=c("educ", "country"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_educ <- ggplot(data=subset(rate_summary_educ, !is.na(educ)), aes(x=educ, y=rate, fill=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("") +
  ylab("") +
  #theme(plot.title = element_text(size = 14, face = "bold"),axis.text.x =element_text(size=14, face ="bold"),axis.text.y =element_text(size=14, face ="bold"),strip.text.x = element_text(size=14, face ="bold"), axis.title.y = element_text(size = 14, face ="bold"))+
  ggtitle("Average rating by education level") +
  scale_fill_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("Low", "Medium", "High")) +
  coord_cartesian(ylim = c(1.27, 6.73)) +
  guides(fill=FALSE)
fig_rate_summary_educ


#km
rate_summary_km <- summarySE(df_conj, measurevar="rate", groupvars=c("country", "km", "id"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_km <- ggplot(data=subset(rate_summary_km, !is.na(km)), aes(x=km, y=rate, color=country, fill=country)) +
  #geom_point() +
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x),
                labels = trans_format("log10", math_format(10^.x))) +
  ggtitle("") +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Yearly driving distance [km]") +
  ylab("Rating") +
  #theme(plot.title = element_text(size = 14, face = "bold"),axis.text.x =element_text(size=14, face ="bold"),axis.text.y =element_text(size=14, face ="bold"),strip.text.x = element_text(size=14, face ="bold"), axis.title.y = element_text(size = 14, face ="bold"))+
  ggtitle("Average rating of policy-packages by yearly driving distance (linear regression)") +
  scale_color_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_fill_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  coord_cartesian(ylim = c(1.27, 6.73)) +
  annotation_logticks(sides = "b")
fig_rate_summary_km


#stopdriving
rate_summary_stopdriving <- summarySE(df_conj, measurevar="rate", groupvars=c("stopdriving", "country"), na.rm = T, conf.interval = 0.95)

fig_rate_summary_stopdriving <- ggplot(data=subset(rate_summary_stopdriving, !is.na(stopdriving)), aes(x=stopdriving, y=rate, color=country, group=country)) + 
  geom_point(stat="identity",
             size=.3) + 
  geom_line() +
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("") +
  ylab("Rating") +
  ggtitle("Average rating of policy-packages by difficulty to stop driving") +
  scale_color_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("Very\ndifficult", "Difficult", "Somewhat\ndifficult", "Neither difficult\nnor easy", "Somewhat\neasy", "Easy", "Very easy", "Non drivers")) +
  coord_cartesian(ylim = c(1.27, 6.73))
fig_rate_summary_stopdriving

##### Anova test for statistical differences ######
library("ggpubr")

#gender
ggboxplot(df_conj, x = c("country", "gender"), y = "rate", 
          color = "country", palette = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"),
          ylab = "rate", xlab = "Groups")

aov_gender <- aov(rate ~ gender + country + gender:country, data = df_conj)
summary(aov_gender)

#educ
ggboxplot(df_conj, x = c("country", "educ"), y = "rate", 
          color = "country", palette = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"),
          ylab = "rate", xlab = "Groups")

aov_educ <- aov(rate ~ educ + country + educ:country, data = df_conj)
summary(aov_educ)

#age
ggboxplot(df_conj, x = c("country", "age"), y = "rate", 
          color = "country", palette = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"),
          ylab = "rate", xlab = "Groups")

aov_age <- aov(rate ~ age + country + age:country, data = df_conj)
summary(aov_age)

#incgrp
ggboxplot(df_conj, x = c("country", "incgrp"), y = "rate", 
          color = "country", palette = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"),
          ylab = "rate", xlab = "Groups")

aov_incgrp <- aov(rate ~ incgrp + country + incgrp:country, data = df_conj)
summary(aov_incgrp)

#joint output
library(export)

filen <- tempfile(pattern = "table_aov")

aov_joint <- aov(rate ~ gender + educ + age + incgrp + country + gender:country + educ:country + age:country + incgrp:country, data = df_conj)

table2spreadsheet(x=NULL,file=filen, type = c("XLS"), digits = 2, digitspvals = 2, add.rownames = TRUE)
summary(aov_joint)
table2spreadsheet(file=filen, width=3.5, font="Calibri", pointsize=12,
          digits=3, digitspvals=1, append=TRUE)

#############paper graphs export #
#conjoint graph all
df_base <- rbind(df_c_base, df_r_base)
df_base$grp <- rep(c("Choice", "Rating"), each=nrow(df_c_base))

p_base <- ggplot(df_base, aes(x=att, y = c, shape = poltyp)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom") +
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) 
p_base


legend4 <- get_legend(p_base)
p_base <- p_base + theme(legend.position="none")

p_base1 <- ggdraw() +
  draw_plot(p_base, x = 0, y = 0.1, width = 1, height = .9) +
  draw_plot(legend4, x = 0, y = 0, width = 1, height = 0.1)
p_base1


ggsave("./figs/p_base_conj_joint1.png", p_base1, width = 9, height = 10)

# Marginal Means graph all
df_mm <- rbind(cmm, rmm)
df_mm$grp <- rep(c("Choice", "Rating"), each=nrow(cmm))

p_mm <- ggplot(df_mm, aes(x=att, y = estimate, shape = poltyp)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  #geom_hline(yintercept = c(0.5, 4), lty="dashed") +
  facet_grid(policy~grp, scales = "free") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull'))
p_mm

legend4 <- get_legend(p_mm)
p_mm <- p_mm + theme(legend.position="none")

p_base1 <- ggdraw() +
  draw_plot(p_mm, x = 0, y = 0.1, width = 1, height = .9) +
  draw_plot(legend4, x = 0, y = 0, width = 1, height = 0.1)
p_base1


ggsave("./figs/marginalmeans_joint.png", p_base1, width = 9, height = 10)

#conjoint graph by country
df_countries <- rbind(df_r_base_US, df_r_base_CN, df_r_base_DE, df_c_base_US, df_c_base_CN, df_c_base_DE)
df_countries$ctr <- rep(c("USA", "China", "Germany", "USA", "China", "Germany"), each=nrow(df_r_base_US))
df_countries$grp <- rep(c("Rating", "Rating", "Rating", "Choice", "Choice", "Choice"), each=nrow(df_r_base_US))
df_countries$poltyp <- factor(ifelse(df_countries$att == "Baseline: No new tax" |  df_countries$att == "Increasing prices by 15%" |  df_countries$att == "Increasing prices by 30%" | 
                                       df_countries$att == "Baseline: General government budget" | df_countries$att == "Public environmental and climate protection programs" | df_countries$att == "Public programs for low-income households" | df_countries$att == "Reduce income taxes" | df_countries$att == "No tax revenues" |
                                       df_countries$att == "Baseline: Keeping subsidies at current level" | df_countries$att == "Halving subsidies" | df_countries$att == "Eliminating subsidies",
                                  "Market-based / push", ifelse(
                                    df_countries$att == "Baseline: Never restricted" | df_countries$att == "1 day per week" | df_countries$att == "3 days per week" | df_countries$att == "5 days per week" | df_countries$att == "Always restricted" | 
                                      df_countries$att == "Baseline: No new requirements" | df_countries$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_countries$att == "Maximum consumption of 4.3 litres per 100 kilometres",
                                    "Non-market-based / push", ifelse(
                                      df_countries$att == "Baseline: No campaigns" | df_countries$att == "Occasional campaigns" | df_countries$att == "Frequent campaigns", 
                                      "Non-market-based / pull", ifelse(
                                        df_countries$att == "Baseline: No increase of support" | df_countries$att == "Reducing prices by 15%" | df_countries$att == "Reducing prices by 30%",
                                        "Market-based / pull", "Other Policy Features")))))


cntrs_output3 <- ggplot(df_countries, aes(x=att, y = c, group = ctr, shape = poltyp, colour = ctr)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom", legend.box = 'vertical') +
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  #ggtitle("Estimated Average Marginal Component Effects") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=c("#e66101", "#5e3c99", "#b2abd2")) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  theme(strip.background = element_rect(fill = 'white'))

cntrs_output2 <- cntrs_output3 + theme(legend.position="none") 

legend <- get_legend(cntrs_output3) 


graph21 <- ggdraw() +
  draw_plot(fig_rate_summary_country, x = 0, y = 0.1, width = .3, height = .25) +
  draw_plot(fig_rate_average, x = .3, y = 0.1, width = .7, height = .25) +
  draw_plot(cntrs_output2, x = 0, y = 0.35, width = 1, height = 0.65) +
  draw_plot(legend, x = 0, y = 0, width = 1, height = 0.1) +
  draw_plot_label(label = c("A", "B", "C"), size = 14,
                  x = c(0, 0, 0.3), y = c(1, 0.35, 0.35))
graph21

ggsave(filename = "./figs/graph21_new.png", graph21, width = 9, height = 13)

#marginal means by country
df_countriesmm <- rbind(rmmUS, rmmCN, rmmDE, cmmUS, cmmCN, cmmDE)
df_countriesmm$ctr <- rep(c("USA", "China", "Germany", "USA", "China", "Germany"), each=nrow(rmmUS))
df_countriesmm$grp <- rep(c("Rating", "Rating", "Rating", "Choice", "Choice", "Choice"), each=nrow(rmmUS))


cntrs_output3 <- ggplot(df_countriesmm, aes(x=att, y = estimate, group = ctr, shape = ctr, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom") +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  #geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=c("#e66101", "#5e3c99", "#b2abd2")) +
  scale_shape_discrete(name="Country") + 
  theme(strip.background = element_rect(fill = 'white'))

cntrs_output2 <- cntrs_output3 + theme(legend.position="none")

legend9 <- get_legend(cntrs_output3)


graph5 <- ggdraw() +
  draw_plot(cntrs_output2, x = 0, y = 0.1, width = 1, height = .9) +
  draw_plot(legend9, x = 0, y = 0, width = 1, height = 0.1) 
graph5

ggsave(filename = "./figs/graphmmCntrs.png", graph5, width = 9, height = 11)

# Marginal Means Choice
df_mm <- rbind(cmm, cmmUS, cmmCN, cmmDE)
df_mm$grp <- rep(c("Pooled", "By country", "By country", "By country"), each=nrow(cmm))
df_mm$ctr <- rep(c("4", "3", "1", "2"), each=nrow(cmm))

choice_mm <- ggplot(df_mm, aes(x=att, y = estimate, group = ctr, shape = poltyp, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom",
        legend.box = "vertical") +
  coord_flip(ylim = c(0.35, 0.65)) +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_y_continuous(breaks = c(0.4, 0.5, 0.6), labels = c("0.4","0.5", "0.6")) +
  scale_color_manual(name="Country", values = c("3" = "#b2abd2", "2" = "#5e3c99", "1" = "#e66101", "4" = "#000000"), labels = c("China", "Germany", "USA", "Pooled")) +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull'))
choice_mm

ggsave("./figs/marginalmeans_choice.png", choice_mm, width = 10.5, height = 10.5)

# Marginal Means Rate
df_mmr <- rbind(rmm, rmmUS, rmmCN, rmmDE)
df_mmr$grp <- rep(c("Pooled", "By country", "By country", "By country"), each=nrow(rmm))
df_mmr$ctr <- rep(c("4", "3", "1", "2"), each=nrow(rmm))

rate_mm <- ggplot(df_mmr, aes(x=att, y = estimate, group = ctr, shape = poltyp, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank(),
        legend.position="bottom",
        legend.box = "vertical") +
  coord_flip(ylim = c(3, 5)) +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_y_continuous(breaks = c(3, 4, 5), labels = c("3","4", "5")) +
  scale_color_manual(name="Country", values = c("3" = "#b2abd2", "2" = "#5e3c99", "1" = "#e66101", "4" = "#000000"), labels = c("China", "Germany", "USA", "Pooled")) +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull'))
rate_mm

ggsave("./figs/marginalmeans_rate.png", rate_mm, width = 10.5, height = 10.5)

#average rating by car dependence
graph3 <- ggdraw() +
  draw_plot(fig_rate_summary_stopdriving, x = 0, y = 0.5, width = 1, height = .5) +
  draw_plot(fig_rate_summary_km, x = 0, y = 0, width = 1, height = .5) +
  draw_plot_label(label = c("A", "B"), size = 14,
                  x = c(0, 0), y = c(1, 0.5))
graph3

ggsave(filename = "./figs/graph31a.png", graph3, width = 8, height = 8)

#Conjoint graph rating including controls (Appendix)
p_r_base_ctrl <- ggplot(df_r_base_ctrl, aes(x=att, y = c, shape = poltyp)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) +
  theme(legend.position="none") +
  facet_grid(policy ~ ., scales = "free")
p_r_base_ctrl

p_r_base_ctrl <- ggdraw() +
  draw_plot(p_r_base_ctrl, x = 0, y = 0.1, width = 1, height = .9) +
  draw_plot(legend4, x = 0, y = 0, width = 1, height = 0.1)
p_r_base_ctrl

ggsave("./figs/p_r_base_ctrl_conj.png", p_r_base_ctrl, width = 9, height = 10)

#average rating by sociodemographics (Appendix)
graphappendix <- ggdraw() +
  draw_plot(fig_rate_summary_gender, x = 0, y = 0.66, width = 0.4, height = .33) +
  draw_plot(fig_rate_summary_educ, x = 0.4, y = 0.66, width = 0.6, height = .33) +
  draw_plot(fig_rate_summary_age2, x = 0, y = 0.33, width = 1, height = .33) +
  draw_plot(fig_rate_summary_incgrp, x = 0, y = 0, width = 0.85, height = .33) +
  draw_plot(legend2, x = 0.85, y = 0.0, width = 0.15, height = 0.33) +
  draw_plot_label(label = c("A", "B", "C", "D"), size = 14,
                  x = c(0, 0.4, 0, 0), y = c(0.99, 0.99, 0.66, 0.33))
graphappendix

ggsave(filename = "./figs/graphappendix.png", graphappendix, width = 9, height = 13.5)


#############paper PPT export
# Marginal Means Choice
df_mm <- rbind(cmm, cmmUS, cmmCN, cmmDE)
df_mm$grp <- rep(c("Pooled", "By country", "By country", "By country"), each=nrow(cmm))
df_mm$ctr <- rep(c("4", "3", "1", "2"), each=nrow(cmm))

choice_mmpp <- ggplot(df_mm, aes(x=att, y = estimate, group = ctr, shape = poltyp, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip(ylim = c(0.35, 0.65)) +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_y_continuous(breaks = c(0.4, 0.5, 0.6), labels = c("0.4","0.5", "0.6")) +
  scale_color_manual(name="Country", values = c("3" = "#b2abd2", "2" = "#5e3c99", "1" = "#e66101", "4" = "#000000"), labels = c("China", "Germany", "USA", "Pooled")) +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull'))
choice_mmpp

ggsave("./figs_retreat/marginalmeans_choicepp.png", choice_mmpp, width = 10, height = 8)

#conjoint graph all
df_base <- rbind(df_c_base, df_r_base)
df_base$grp <- rep(c("Choice", "Rating"), each=nrow(df_c_base))

p_base <- ggplot(df_base, aes(x=att, y = c, shape = poltyp)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_shape_manual(name="Policy type", values = c(18:15), breaks = c('Market-based / push','Non-market-based / push','Market-based / pull','Non-market-based / pull')) 
p_base

ggsave("./figs_retreat/retreat_fig1.png", p_base, width = 10, height = 8)

#AMCE by country
df_countries <- rbind(df_r_base_US, df_r_base_CN, df_r_base_DE, df_c_base_US, df_c_base_CN, df_c_base_DE)
df_countries$ctr <- rep(c("USA", "China", "Germany", "USA", "China", "Germany"), each=nrow(df_r_base_US))
df_countries$grp <- rep(c("Rating", "Rating", "Rating", "Choice", "Choice", "Choice"), each=nrow(df_r_base_US))


fig3_country <- ggplot(df_countries, aes(x=att, y = c, group = ctr, shape = ctr, colour = ctr)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  #ggtitle("Estimated Average Marginal Component Effects") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=c("#e66101", "#5e3c99", "#b2abd2")) +
  scale_shape_discrete(name="Country") + 
  theme(strip.background = element_rect(fill = 'white'))

ggsave(filename = "./figs_retreat/fig3_country.png", fig3_country, width = 10, height = 8)

#marginal means by country
df_countriesmm <- rbind(cmmUS, cmmCN, cmmDE)
df_countriesmm$ctr <- rep(c("USA", "China", "Germany"), each=nrow(cmmUS))
df_countriesmm$grp <- rep(c("Choice", "Choice", "Choice"), each=nrow(cmmUS))


cntrs_output3 <- ggplot(df_countriesmm, aes(x=att, y = estimate, group = ctr, shape = ctr, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 0.5, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=c("#e66101", "#5e3c99", "#b2abd2")) +
  scale_shape_discrete(name="Country") + 
  theme(strip.background = element_rect(fill = 'white'))


ggsave(filename = "./figs_retreat/graphmmCntrs.png", cntrs_output3, width = 8, height = 8)

#average rating
rateavrg1 <- aggregate(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                        attrib5_lab +  attrib6_lab +  attrib7_lab, df_conj, mean)

fig_rate_averageall <- ggplot(rateavrg1, aes(x=rate)) + 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" ),
        axis.text=element_text(size=12), axis.title=element_text(size=16)) +
  xlab("Distribution of average ratings for policy-packages") +
  ylab("Density") +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  coord_cartesian(ylim = c(.03095, 0.65), xlim = c(1.27, 6.73)) +
  guides(colour=FALSE) +
  geom_vline(xintercept = 4, lty="dashed")

ggsave(filename = "./figs_retreat/avrgratingall.png", fig_rate_averageall, width = 7, height = 6)

fig_rate_average1 <- ggdraw() +
  draw_plot(fig_rate_average, x = 0, y = 0, width = 0.85, height = 1) +
  draw_plot(legend3, x = 0.9, y = 0, width = 0.1, height = 1)
fig_rate_average1

ggsave(filename = "./figs_retreat/avrgrating.png", fig_rate_average1, width = 8, height = 6)


fig_rate_average2 <- ggplot(rateavrg, aes(x=rate, colour=country)) + 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80"),
        legend.text=element_text(size=12),
        axis.text=element_text(size=12), axis.title=element_text(size=16)) +
  xlab("Distribution of average ratings for policy-packages") +
  ylab("Density") +
  scale_colour_manual(name="", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  coord_cartesian(ylim = c(.0257, 0.55), xlim = c(1.27, 6.73)) +
  geom_vline(xintercept = 4, lty="dashed")


ggsave(filename = "./figs_retreat/avrgratingctrs.png", fig_rate_average2, width = 8, height = 6)

#marginal means Rating
df_countriesmm <- rbind(rmmUS, rmmCN, rmmDE, rmm)
df_countriesmm$ctr <- rep(c("USA", "China", "Germany", "Pooled"), each=nrow(cmmUS))
df_countriesmm$grp <- rep(c("Rating", "Rating", "Rating", "Rating"), each=nrow(cmmUS))

cols <- c("China" = "#e66101","Germany" = "#5e3c99","USA" = "#b2abd2", "Pooled" = "#000000")
shap <- c("China" = 16,"Germany" = 17,"USA" = 15, "Pooled" = 3)

mmrating <- ggplot(df_countriesmm, aes(x=att, y = estimate, group = ctr, shape = ctr, colour = ctr)) +
  geom_pointrange(aes(min = estimate - 1.95 * std.error, max = estimate + 1.95 * std.error), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip() +
  ylab("Marginal Means") +
  xlab("Attributes") +
  geom_hline(yintercept = 4, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=cols, breaks = c('China','Germany','USA','Pooled')) +
  scale_shape_manual(name="Country", values=shap, breaks = c('China','Germany','USA','Pooled')) + 
  theme(strip.background = element_rect(fill = 'white'))


ggsave(filename = "./figs_retreat/mmrating.png", mmrating, width = 8, height = 8)
